Skip to content

Defect repairs #1432

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Merged
jeffriley merged 4 commits into from
Aug 23, 2025
Merged

Defect repairs #1432

jeffriley merged 4 commits into from
Aug 23, 2025

Conversation

@jeffriley
Copy link
Collaborator

@jeffriley jeffriley commented Aug 21, 2025

@ilyamandel a couple more minor fixes related to Hurley coefficients and constants:

  • Clamped B_GAMMA to [0.0, B_GAMMA] (per discussion just after eq 23 - confirmed in BSE Fortran source)
  • Corrected calculation for Hurley Gamma constant C (C_GAMMA - see Hurley et al. 2000, just after eq 23, should use a(75) <= 1.0, not a(75) == 1.0 - confirmed in BSE Fortran source)

@jeffriley jeffriley requested a review from ilyamandel August 21, 2025 05:48
@jeffriley jeffriley added bug Something isn't working severity_minor This bug is not very severe urgency_low This issue is not urgent labels Aug 21, 2025
ilyamandel
ilyamandel approved these changes Aug 21, 2025
View reviewed changes
Copy link
Collaborator

@ilyamandel ilyamandel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can I suggest that we reference the Fortran code in the comments along with Hurley Eq. (23) -- e.g., "// see Hurley et al. 2000, eq 23, but the code here follows the Fortran SSE implementation". Because the previous version implemented Eq. 23 as written in the paper, the new one does not (but I trust you that it implements the Fortran code).

@jeffriley
Copy link
Collaborator Author

jeffriley commented Aug 21, 2025 •
edited
Loading

Because the previous version implemented Eq. 23 as written in the paper, the new one does not

Do you mean the addition of the max()? That's in the discussion following eq 23 (as is the <= 1.0)

Edit: When I first made the change I added a max() to C_GAMMA (to ensure it is >= 0.0), but I wasn't sure the comment didn't refer only to B = gamma (when M = 1.0)... Checking MainSequence::CalculateGamma(const double p_Mass) now, we don't ensure gamma is >= 0.0 - I'm pretty sure Hurley does in the BSE code - I'll confirm that tomorrow.

@jeffriley
Copy link
Collaborator Author

jeffriley commented Aug 21, 2025 •
edited
Loading

Here's the state of affairs:

Current COMPAS code: Mainsequence.cpp::CalculateGamma(), BaseStar::CalculateAnCoefficients() (GammaConstants(B_GAMMA) and GammaConstants(C_GAMMA))
Legacy code: star.cpp, calculate_gamma(), calculate_B_gamma(), calculate_C_gamma()

BSE code zfuncs.f:

*
* A function to evaluate the radius gamma coefficent.
* (JH 24/11/97)
*
m1 = 1.d0
if(m.gt.(a(88)+0.1d0))then                                            <-- see discussion immediately prior to eq 23; endpoint is 0.1, not 1.0
    rgammf = 0.d0                                                     <-- gamma = 0.0
else
    b1 = MAX(0.d0,a(83) + a(84)*(m1-a(85))**a(86))                    <-- this is B (in paper) (B_GAMMA in COMPAS) (but is also used as C (in paper, C_GAMMA in COMPAS)).  MAX() missing in legacy code and current COMPAS code
    if(m.le.m1)then                                                   <-- eq 23 condition M <= 1.0
        rgammf = a(83) + a(84)*ABS(m-a(85))**a(86)                    <-- gamma - ABS() missing in legacy code and current COMPAS code
    elseif(m.le.a(88))then                                            <-- eq 23 condition 1.0 < M <= a[75]
        rgammf = b1 + (a(89) - b1)*((m - m1)/(a(88) - m1))**a(87)     <-- gamma
    else                                                              <-- eq 23 condition a[75] < M < a[75] + 0.1 (note 0.1, not 1.0)
        if(a(88).gt.m1) b1 = a(89)                                    <-- C = a[80] (C_GAMMA in COMPAS) - see discussion following eq 23; note reuse of b1 as C (was B)
        rgammf = b1 - 10.d0*b1*(m - a(88))                            <-- gamma
    endif
    rgammf = MAX(rgammf,0.d0)                                         <-- ensure gamma +ve.  MAX() missing in current COMPAS code; exists in legacy code
endif

return
end
    

BSE code a(83) = paper a[76]
BSE code a(84) = paper a[77]
BSE code a(85) = paper a[78]
BSE code a(86) = paper a[79]
BSE code a(87) = paper a[81]
BSE code a(88) = paper a[75]
BSE code a(89) = paper a[80]

So, I will push another change that adds ABS() appropriately, and clamps gamma to [0.0, gamma], in Mainsequence.cpp::CalculateGamma()

ilyamandel
ilyamandel reviewed Aug 23, 2025
View reviewed changes
Copy link
Collaborator

@ilyamandel ilyamandel left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

GammaConstants(C_GAMMA) = (utils::Compare(a[75], 1.0) <= 0) ? GammaConstants(B_GAMMA) : a[80]; // Hurley et al. 2000, eq 23 and discussion immediately following - <= 1.0 confirmed in BSE Fortran code

But the Hurley paper actually says
"and C=a80 unless a75=1.0 when C=B"
Equal, not \leq.
That's why I was suggesting referencing the Fortran code directly in opposition to the paper -- but I am OK with the current comment wording if you prefer that.

@jeffriley
Copy link
Collaborator Author

jeffriley commented Aug 23, 2025

"and C=a80 unless a75=1.0 when C=B"

Not the version I'm looking at...

hurleyp9.pdf

@jeffriley jeffriley merged commit d224356 into dev Aug 23, 2025
3 checks passed
@jeffriley jeffriley deleted the minor-fix branch August 23, 2025 00:23
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@ilyamandel ilyamandel ilyamandel approved these changes

Assignees

No one assigned

Labels

bug Something isn't working severity_minor This bug is not very severe urgency_low This issue is not urgent

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

3 participants

@jeffriley @ilyamandel